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K^ ■ Abstract 

j^ ' Galaxies are arranged in interconnected walls and filaments forming a cosmic web en- 

JH I compassing huge, nearly empty, regions between the structures. Many statistical methods 

■ - - ' have been proposed in the past in order to describe the galaxy distribution and discriminate 

the different cosmological models. We present in this paper results relative to the use of 
new statistical tools using the 3D isotropic undecimated wavelet transform, the 3D ridgelet 
transform and the 3D beamlet transform. We show that such multiscale methods produce 
a new way to measure in a coherent and statistically reliable way the degree of clustering, 
filamentarity, sheetedness, and voidedness of a dataset. 
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1 Introduction 

Galaxies are not uniformly distributed throughout the universe. Voids, filaments, clusters, and 
walls of galaxies can be observed, and their distribution constraints our cosmological theories. 



Therefore we need reliable statistical methods to compare the observed galaxy distribution with 
theoretical models and cosmological simulations. 

The standard approach for testing models is to define a point process which can be char- 
acterized by statistical descriptors. This could be the distribution of galaxies of a specific type 
in deep redshift surveys of galaxies (or of clusters of galaxies). In order to compare models of 
structure formation, the different distribution of dark matter particles in N-body simulations 
could be analyzed as well, with the same statistics. 

The two-point correlation function ^(r) has been the primary tool for quantifying large-scale 
cosmic structure |24| . Assuming that the galaxy distribution in the Universe is a realization 
of a stationary and isotropic random process, the two-point correlation function can be defined 
from the probability 5P of finding an object within a volume element 6V at distance r from 
a randomly chosen object or position inside the volume: 6P = n(l + i{r))5V ^ where n is the 
mean density of objects. The function ^(r) measures the clustering properties of objects in a 
given volume. It is zero for a uniform random distribution, positive (respectively, negative) for 
a more (respectively, less) clustered distribution. For a hierarchical clustering or fractal process, 
1 + ^(r) follows a power-law behavior with exponent D2 — 3. Since ^(r) ~ r~^ for the observed 
galaxy distribution, the correlation dimension for the range where ^(r) S> 1 is D2 — 3 — 7. The 
Fourier transform of the correlation function is the power spectrum. The direct measurement 
of the power spectrum from redshift surveys is of major interest because model predictions are 
made in terms of the power spectral density. It seems clear that the real space power spectrum 
departs from a single power-law ruling out simple unbounded fractal models [HEl- The two- 
point correlation function can been generalized to the N-point correlation function j35l I25j . 
and all the hierarchy can be related with the physics responsible for the clustering of matter. 
Nevertheless they are difficult to measure, and therefore other related statistical measures have 
been introduced as a complement in the statistical description of the spatial distribution of 
galaxies J2D|, such as the void probability function |2Tj, the multifractal approach ^H]; the 
minimal spanning tree ^ 1151 [S] , the Minkowski functionals j22l I13j or the J function |171 I14j 
which is defined as the ratio J{r) = ^prl , where F is the distribution function of the distance 

r of an arbitrary point in R to the nearest object in the catalog, and G is the distribution 
function of the distance r of an object to the nearest object. Wavelets have also been used for 
analyzing the projected 2D or the 3D galaxy distribution [HI OHIIT^ OHIIT^ . 

The Genus statistic flO] measures the topology or the degree of connectedness of the un- 
derlying density field. Phase correlations of the density field can be detected by means of this 
topological analysis. The Genus is calculated by (i) convolving the data by a kernel, generally a 
Gaussian, (ii) setting to zero all values under a threshold v in the obtained distribution, and (iii) 
taking the difference D between the number of holes and the number of isolated regions. The 
Genus curve G{v) is obtained by varying the threshold level v. The first step of the algorithm, 
the convolution by a Gaussian, may be dramatic for the description of filaments, which are 
spread out along all directions, as it will be shown in next section. The Genus is related with 
one of the four Minkowski functionals that describe well the overall morphology of the galaxy 
distribution ^^ . Minkowski functionals have been used to elaborate sophisticated tools to mea- 
sure the filamentarity and planarity of the distribution by means of shape finders quantities 

m- 

The Sloan Digital Sky Survey (Early Data Release) has recently be analyzed using a 3D 
Genus Statistics ^J and results were consistent with that predicted by simulations of a A- 
dominated spatially-flat cold dark matter model. 

New multiscale methods have recently emerged, the beamlet transform jUEj and the ridgelet 
transform [3], which allows us to better represent data containing respectively filaments and 
sheets, while wavelets represent well isotropic features (i.e. cluster in 3D). As each of these 



three transforms represents perfectly one kind of feature, all of them are useful and should be 
used to describe a given catalog. 

Section 2 describes the 3D wavelet transform, and how wavelets can be used for estimating 
the underlying density. Sections 3 and 4 describes respectively the 3D ridgelet transform and the 
3D beamlet transform. It is shown in section 4 through a set of of experiments how these three 
3D transforms can be combined in order to describe statistically the distribution of galaxies. 

2 The 3D Wavelet Transform 

2.1 The Undecimated Isotropic Wavelet Transform 

The wavelet transform of a signal produces, at each scale j, a set of zero-mean coefficient values 
{wj}. Using an algorithm such as the undecimated isotropic wavelet decomposition [SSI, ^^i^ set 
{wj} has the same number of pixels as the signal and thus this wavelet transform is a redundant 
one. Furthermore, using a wavelet defined as the difference between the scaling functions of two 
successive scales 

the original cube c = cq can be expressed as the sum of all the wavelet scales and the smoothed 
array cj 

J 

^0,x,y,z — Cj^x,y,z ~r / ^ Wj^x,y,z ; \^) 

i=i 

The set w = {wi , W2, ...,wj,cj}, where cj is a last smooth array, represents the wavelet transform 
of the data. If we denote as W the wavelet transform operator and A^ the pixel number of c, the 
wavelet transform w {w = Wc) has (J + 1)A'^ pixels (redundancy factor of J + 1). The scaling 
function (f) is generally chosen as a spline of degree 3, and the 3D implementation is based on 
three ID sets of (separable) convolutions. Like the scaling function (/>, the wavelet function -0 is 
isotropic (point symmetric). More details can be found in [SSI EH- 
For each a > 0, 6i, 621 &3 £ R'^ 1 the wavelet is defined by 

i'aMMM '■ R ^ R 

Given a function / € L^(R'^), we define its wavelet coefficients by: 
W/ : R^ ^ R _ 

VVfia, 61, 62, ^3) = / 0a,fcj,fc2,b3(x)/(x)dx. 
Figure n shows an example of 3D wavelet function. 

2.2 3D Galaxy Distribution Filtering 

For the noise model, given that this relates to point pattern clustering, we have to consider the 
Poisson noise case. The autoconvolution histogram method used for X-ray image |34j can also 
be used here. It consists to calculate numerically the probability distribution function (pdf) 
of a wavelet vjj^x,y,z coefficient with the hypothesis that the galaxies used for obtaining Wj^x,y,z 
are randomly distributed. The pdf is obtained by autoconvolving n times the histogram of the 
wavelet function, n being the number of galaxies which have been used for obtaining Wj^x,y,z, 
i.e. the number of galaxies in a box around (x, y, x), the size of the box depending of the scale 
j. More details can be found in jHHI^ . 
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Figure 1: Example of wavelet function. 



Once the pdf relative to the coefficient Wj^x,y,z is known, we can detect the significant wavelet 



coefficients easily. We derive two threshold values T, 



--,, and r,--,, such that 



Prob{W<T^:i,) = e 

Prob{W>T;^:i,) = e (3) 

e corresponding to the confidence level, and the positive (respective negative) wavelet coefficient 
is significant if it is larger than T"!"^ ^ (resp. lower than TJ™" ,). 

A simple filtering method would now consist to set to zero (i.e. thresholding) all unsignificant 
coefficients and reconstruct the filtered data cube by addition of the different scales. But when 
a redundant wavelet transform is used, the result after a simple hard thresholding can still be 
improved by iterating [50]. We want the wavelet transform of our solution s to reproduce the 
same significant wavelet coefficients (i.e., coefficients larger than Tj). This can be expressed in 
the following way: 



{Ws)j^k = Wj^k if Wj^k is significant 



(4) 



where Wj^k ^-re the wavelet coefficients of the input data s at scale j and at position k = (x, y, z). 
The relation is not necessarily verified in the case of non-orthogonal transforms, and the resulting 



effect is generally a loss of flux inside the objects. The residual signal (i.e. s — s) still contains 
some information at positions where the objects are. 

Denoting M the multiresolution support of s (i.e. Mj^ = 1 if Wj^k is significant, and 
otherwise), we want: 

M.Ws = M.Ws 

The solution can be obtained by the following Van Cittert iteration j33j : 

gu+i ^ s" + W-^M.Ws - M.Ws") 

= s" + W-^M.Wi?") (5) 

where K^ = s — s". 

Iterative Filtering with a Smoothness Constraint 

A smoothness constraint can be imposed on the solution. 

min||>Vs||£^, subject to s £ C, (6) 

where C is the set of vectors s which obey the linear constraints 

Sk > 0,\/k , . 

I (Ws-W%fc|<e,-; ^'^ 

Here, the second inequality constraint only concerns the set of significant coefficients, i.e. those 
indices which exceed (in absolute value) a detection threshold tj. Given a tolerance vector 
e = {ei, ...,ej}, we seek a solution whose coefficients (Ws)j_fc, at scale and position where 
significant coefficients were detected, are within Cj of the noisy coefficients {Ws)j^k- For example, 
we can choose Cj = (Tj/2. In short, our constraint guarantees that the reconstruction be smooth 
but will take into account any pattern which is detected as significant by the wavelet transform. 
We use an £i penalty on the coefficient sequence because we are interested in low complex- 
ity reconstructions. There are other possible choices of complexity penalties; for instance, an 
alternative to ^ would be 

min||s||Ty, subject to s £ C. 

where || • \\tv is the Total Variation norm, i.e. the discrete equivalent of the integral of the 
Euclidean norm of the gradient. Expression © can be solved using the method of hybrid 
steepest descent (HSD) fSH'. HSD consists of building the sequence 



~n+l 



P(s") - A„+iVj(P(5")); (8) 



here, P is the £2 projection operator onto the feasible set C, Vj is the gradient of equation ©, 
and (A„)n>i is a sequence obeying (An)n>i G [0, 1] and lim„^+oo A^ = 0. 

Unfortunately, the projection operator P is not easily determined and in practice we use the 
following proxy; compute Ws and replace those coefficients which do not obey the constraints 
I (Ws — \Vs)j^k\ < Cj (those which fall outside of the prescribed interval) by those of s; apply the 
inverse transform. 

The filtering algorithm is: 



1. Liitialize Lmax = 1; the number of iterations Ni, and 5> 



A^, 



2. Estimate the noise standard deviation o"s, and set Cj = -^ for all j. 

3. Calculate the transform: w^'^' = Ws. 



4. Set X = Lr, 



n 



0, and s" to 0. 



5. While A >= do 



• u 



Calculate the transform a = Wtt. 
For all coefficients Uj^k do 



* Calculate the residual r 



'j,k 



W 



j,k 



O^i.k 



""Ik 



* if w- ^ is significant and | rjj. \ > Cj then Uj^k 

* ai,fc = sgn{aj^k){\ ctj^k I -Ao-j)+- 
- u = W^^a 

• Threshold negative values in u and s""'"^ = u. 

• n = n + 1, A = A — (5a, and goto 5. 

In practice, a small number of iterations (<10) is enough. The (.)+ operator means that negative 
values are set to zero ((a)+ = MAX(0, a)). 

Experiments 
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Figure 2: Simulated data. 

Fig. [21 shows a simulation of a 60/i~^ Mpc box of universe, made by A. Klypin (simulated 
data at http://astro.nmsu.edu/~aklypin/PM/pmcode). It represents the distribution of dark 
matter in the present-day universe and each point is a dark matter halo where visible galaxies 



are expected to be located, i.e. the distribution of dark matter halos can be compared with the 
distribution of galaxies in catalogs of galaxies. 






Figure 3: Simulated data filtered by a Gaussian filter (two left panels corresponding respectively 
to o" = 1 and a = 3) and using the wavelet algorithm described in the text (right panel). 
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Figure 4: The genus curve for the N-body model shown in Fig. El convolving the data with a 
Gaussian filter with different values of a and the genus for the wavelet filtered data set. 

Fig.Elshows, at the right panel, the same data set filtered by the 3D wavelet transform, using 
the algorithm described previously. The two left panels correspond to Gaussian smoothing 



Ty(x) 



1 



(27r)3/V3 



exp 



X 



2a2)' 



(9) 



with (7 = 1 and a = 3. We can see that when the bandwidth is too small the discreteness and 
the noise dominate the density reconstructed field, while large value of a erase all the small 
scale features of the distribution. This is illustrated in Fig. \^ where we can see the strong 
dependence of the genus curve with the width of the Gaussian filter, being its interpretation 
completely dependent of the choice of the bandwidth. The wavelet reconstructed density field 
keeps information at all scales due to its multiscale nature. This is observed in the 3D image 
at the right panel of Fig. El where we see how large filaments, big clusters and walls coexist 



with small scale features such as the density enhancement around groups and small clusters. 
The genus curve of this adaptive reconstructed density field is much more informative because 
it does not depend of the particular choice of the filter radius. 

3 The 3D Ridgelet Transform 

3.1 The 2D Ridgelet Transform 

The two-dimensional continuous ridgelet transform of a function / G L^(R^) is defined as follows 

We first select a smooth function ij) G L^(R), we assume that ijj satisfies the admissibility 
condition 



mowmdc <oc, 



(10) 



which holds if -0 has a sufficient decay and a vanishing mean j ip{t)dt = Q {ip can be normalized 
so that it has unit energy l/(27r) J" |V'(?)P'^'? = !)• 
For each a > 0, 6 G R and 9i G [0, 27r[, we define the ridgelet by 

V'a,fe,6»i : R ^ R 

V'a,fe,6»i(a;i,X2) = a"^/^ • il){{xicos6i + X2sin6'i - 6)/a); 
Given a function / G L^(R^), we define its ridgelet coefficients by: 

7^/ : R3 ^ R _ 

%(a, b,ei) = J V'a,b,ei (x)/(x)dx. 
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Figure 5: Definition of anglel 0\ and Q2 in R^ and R^ 

It has been shown |;3^ that the ridgelet transform is precisely the application of a 1-dimensional 
wavelet transform to the slices of the Radon transform (where the angular variable Q\ is con- 
stant). This method is therefore optimal to detect lines of the size of the image (the integration 
increase as the length of the line). More details on the implementation of the digital ridgelet 
transform can be found in j31j . 

Figure IHl (left) shows an example ridgelet function. This function is constant along lines 
x\ cos d + X2 sin = const. Transverse to these ridges it is a wavelet (see figure El (right) . 



3.2 From 2D to 3D 

The three-dimensional continuous ridgelet transform of a function / G L^(R^) is given by: 




Figure 6: Example of 2D ridgelet function. 

7^/ : R^ ^ R 

TZfia, b, 6*1, 6*2) = / V'a,6,ei,e2 W/W^^- 
where a > 0, 6 G R , 6*1 G [0, 27r[ and 02 G [0, tt[. 
The ridgelet function is defined by: 
'ipa,b,ei,e2 ■ R ^ R 

i^a,b,ei,92{xi^X2,X3) = a^^l"^ ■ V'((xicos6'icos6'2 + X2 sin 6*1 cos 6*2 + 2:3 sin 6*2 - ft)/a); 
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Figure 7: Example of ridgelet function. 

Figure shows an example of ridgelet function. It is a wavelet function in the direction 
defined by the line {9i,theta2), and it is constant along the orthogonal plane to this line. 

As in the 2D case, the 3D ridgelet transform can be built by extracting lines in the Fourier 
domain. Let 0(11,12,13) a cube of size {N,N,N), the algorithm consists in the following steps: 



1. 3D-FFT. Compute c{ki, k2, k^), the three-dimensional EFT of the cube c(ii, ^2; ^3)- 
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Figure 8: 3D ridgelet transform flowgraph. 

2. Cartesian to Spherical Conversion. Using an interpolation scheme, substitute the sampled 
values of c obtained on the Cartesian coordinate system (/ci,fc2,/c3) with sampled values 
of on a spherical coordinate system {9i, O2, p)- 

3. Extract lines. Extract the 3A^^ lines (size A^) passing through the origin and the boundary 
of c. 

4. ID-IFFT. Compute the one-dimensional inverse FFT on each line. 

5. ID-WT. Compute the one-dimensional wavelet transform on each line. 

Figure |H1 the 3D ridgelet transform flowgraph. The 3D ridgelet transform allows us to detect 
sheets in a cube. 

Local 3D Ridgelet Transform 

The ridgelet transform is optimal to find sheets of the size of the cube. To detect smaller 
sheets, a partitioning must be introduced [2]. The cube c is decomposed into blocks of lower 
side-length h so that for a N * N * N cube, we count N/b blocks in each direction. After the 
block partitioning. The detection is therefore optimal for sheets of size 6x6 and of thickness aj , 
Oj corresponding to the different dyadic scales used in the transformation. 

4 The 3D Beamlet Transform 

4.1 Definition 

The X-ray transform of a continuum function f{x, y, z) with (x, y, z) G R'^ is defined by 

{Xf){L) = j f{p)dp (11) 

where L is a line in R^, and p is a variable indexing points in the line. The transformation 
contains all line integrals of /. The Beamlet Transform (BT) can be seen as a multiscale digital 
X-ray transform. It is multiscale transform because, in addition to the multiorientation and 
multilocation line integral calculation, it integrated also over line segments at different length. 
The 3D BT is an extension to the 2D BT, proposed by Donoho and Huo |4j. 
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The system of 3D beams 

The first choice to consider is the hne segment set. We would like to have an expressive set 
of line segments in the sense that it includes line segments with various lengths, locations and 
orientations lying inside a 3D volume and at the same time has reasonable size. 

A seemingly natural candidate for the set of line segments is the family of all line segments 
between any voxel corner and any other voxel corner, the set of 3-D beams. The beams set 
is expressive but can be of huge cardinality for even moderate resolutions. For a 3D data set 
with rfi voxels we get 0{n^) 3D beams - So that is clearly infeasible to use the collection of 3-D 
beams as a basic data structure since any algorithm based on this set will have a complexity 
with lower bound of n^ and hence unworkable for typical sizes 3-D images. 



4.2 The Beamlet System 

A dyadic cube C(/ci, A;2, ^3, j) C [0, 1]^ is the collection of points 

{(xi,X2,X3) : [ki/2^, {ki + l)/2^'] X [fc2/2^ (fca + l)/2^'] x [fc3/2^ {k^ + l)/2^']} 

where < fci , A;2 , A;3 < 2-' for an integer j > 0. We will refer to j as the scale of the dyadic cube 
Such cubes can be viewed as descended from the unit cube C(0, 0, 0, 0) = [0, l]'^ by recursive 
partitioning. Hence, the result of splitting C(0, 0, 0, 0) in half along each axis is the eight cubes 
C(A;i, A;2, fes, 1) where ki G {0, 1}, splitting those in half along each axis we get the 64 subcubes 
C{ki,k2,k3,2) where ki G {0,1,2,3}, and if we decompose the unit cube into n^ voxels using 
a uniform n-by-n-by-n grid with n = 2'^ dyadic, then the individual voxels are the n^ cells 
C{ki, k2, ks, J), < ki, k2, k^ < n. 
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Figure 9: Dyadic cubes 

Associated to each dyadic cube we can build a system of line segments that have both of their 
end-points lying on the cube boundary. We call each such segment A beamlet. If we consider 
all pairs of boundary voxel corners we get 0{n'^) beamlets for a dyadic cube with a side length 
of n voxels, we will work with a slightly different system in which each line is associated with 
a slope and an intercept instead of its end-points as will be explained below. However, we will 
still have O(n^) cardinality. Assuming a voxel size of 1/n we get J + 1 scales of dyadic cubes 
where n = 2'^ , for any scale < j < J there are 2^^ dyadic cubes of scale j and since each dyadic 
cube at scale j has a side length of 2 -^ voxels we get 0(2^*^ •'•*) beamlets associated with the 
dyadic cube and a total of 0(2^'^^^) = 0{n^/2^) beamlets at scale j. If we sum the number of 
beamlets at all scales we get 0{n'^) beamlets. 

We have constructed above a multi-scale arrangement of line segments in 3D with controlled 
cardinality of 0(71"^), the scale of a beamlet is defined as the scale of the dyadic cube it belongs 
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to so lower scales correspond to longer line segments and finer scales correspond to shorter line 
segments. Figure ITUl shows 2 beamlets at different scales. 





Figure 10: Examples of Beamlets at two different scales, (a) Scale (coarsest scale) (b) Scale 1 
(next finer scale). 

To construct the set of beamlets for a given dyadic cube we use the slope-intercept pairs 
system. For a data cube of n x n x re voxels consider a coordinate system with the cube center 
of mass at the origin and a unit length for a voxel. Hence, for (x, y, z) in the data cube we 
have \x\, \y\, \z\ < n/2. We can consider three kinds of lines: x-driven, y-driven, and z-driven, 
depending on which axis provides the shallowest slopes. An j;-driven line takes the form 



SzX + t. 



y 



SyX ~r f^y 



with slopes Sz,Sy, and intercepts tz and ty. Here the slopes \sz\, \sy\ < 1. y- and z-driven lines 
are defined with an interchange of roles between x and y or z, as the case may be. 

We will consider the family of lines generated by this, where the slopes and intercepts run 
through an equispaced family: 

Sx, Sy, Sz G {2£/re : I = -re/2, . . . , n/2 - 1}, tx,ty,tz G {^ : -re/2, . . . , n/2 - 1}. 

Using the above family of lines for a given data cube we can define the set of beamlets belong 
to the cube to be the set of line segment obtained by taking the intersection of each line with 
the cube. 

With the choice of indices range above we get that all beamlets associated with a data cube 
of size re have lengths bigger than re/2, half of the cube length and \/3n, the cube main diagonal 
length. 

Computational aspects 

The Beamlet coefficients are the line integrals over the set of beamlets. A digital 3-D image can 
be regarded as a 3-D piece-wise constant function and each line integral is just a weighted sum 
of the voxel intensities along the corresponding line segment. Donoho and Levi ^0. discuss in 
detail different approaches for computing line integrals in a 3-D digital image. Computing the 
beamlet coefficients for real applications data sets can be a challenging computational task since 
for a data cube with n x n x n voxels we have to compute O(re^) coefficients. By developing 
efficient cache aware algorithms we are able to handle 3-D data sets of size up to n = 256 on a 
single fast machine in less than a day running time. We will mention that in many cases there 
is no interest in the coarsest scales coefficient that consumes most of the computation time and 
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in that case the over ah running time can be significantly faster. The algorithms can also be 
easily implemented on a parallel machine of a computer cluster using a system such as MPI in 
order to solve bigger problems. 

4.3 The FFT-based transformation 

Let ip S -L^(R^) a smooth function satisfying the admissibility condition (a 2D wavelet function), 
the three-dimensional continuous beamlet transform of a function / G L^(R^)is given by: 
Bf. R5 ^ R 

Bf{a,bi,b2, 01,62) = /V'a,6,0i,e2W/W^^- 
where a > 0, 61,62 G R ,6*1 G [0,27r[ and 02 G [0,7r[. 
The beamlet function is defined by: 
'ipa,bi,b2,ei,e2 '■ R ^ R 



i^a,hub2,0i,e2{xi,X2,X3) 



-1/2. 



ip{ {—xism6i + X2Cos9i + bi)/a, 

{xi cos 6*1 cos 02 + X2 sin 01 cos 02 — 2:3 sin 02 + b2)/a); 
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Figure 11: Example of beamlet function. 

Figured shows an example of beamlet function. It is constant along lines of direction (0i, 
92), and a 2D wavelet function along plane orthogonal to this direction. 
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Figure 12: 3D beamlet transform flowgraph. 

The 3D beamlet transform can be built using the "Generalized projection-slice theorem" 
[Sni. Let 

• /(x) an n dimensional function, 

• TZadmf its m-dimensional partial radon transform along the first m cardinal directions, 
m < n, TZadmf is a function of (p, /i^; Xm+i, ■■■,Xn), fJ-m ^ unit directional vector in TZadm 
(note that for a given projection angle, the m dimensional partial radon transform of /(x) 
has (n — m) untransformated spatial dimension and a (n-m+l) dimensional projection 
profile) , 

• {J^/}(k) its Fourier transform (x and k are conjugate variable pairs of ^), 

The Fourier transform of the m dimensional partial radon transform IZmf is related to the 
Fourier transform of / {J-f) by the following relation 

{^n-m+lT^admf}{k, km+1, ■■■, K) = {J^f}{k^l^, km+1, ■■■, kn) (12) 

Let c{ii,i2,i'i) a cube of size {N, N, N), the Beamlet algorithm consists in the following steps: 

1. 3D-FFT. Compute c{ki, /c2, k^), the three-dimensional FFT of the cube c(ii, ^2, ^3)- 

2. Cartesian to Spherical Conversion. Using an interpolation scheme, substitute the sampled 
values of c obtained on the Cartesian coordinate system [ki,k2,k^) with sampled values 
of on a spherical coordinate system {9i, 62, p)- 

3. Extract planes. Extract the 3A^^ planes (of size A^ x A^) passing through the origin (each 
line used in the 3D ridgelet transform defines set of orthogonal planes, we take the one 
which include the origin). 

4. 2D-IFFT. Compute the two-dimensional inverse FFT on each plane. 

5. 2D-WT. Compute the two-dimensional wavelet transform on each plane. 

Figure [T2l the 3D beamlet transform flowgraph. The 3D beamlet transform allows us to detect 
filament in a cube. The beamlet transform algorithm presented in this section differs from the 
one presented in Q, and relation between both algorithms is given in 0. 
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5 Experiments 

5.1 Experiment 1 

We have simulated three data set containing respectively a cluster, a plane and a line. On each 
data set, Poisson noise have been added with eight different background levels. Then we have 
applied the three transforms on the 24 simulated data set. The coefficients distribution related 
to each transformation is normalized using twenty realizations of a 3D flat distribution with a 
Poisson noise and which have the same number of counts as in the data. 
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Figure 13: Simulation of cubes containing a cluster (top), a plane (middle) and a line (bottom). 

Figure ^J shows, from top to bottom, the maximum value of the normalized distribution 
versus the noise level for our three simulated data set. As expected, wavelets, ridgelets and 
beamlets are respectively the best for clusters, sheets and lines detection. It must also be 
underlined that a feature can be detected with a very high signal-to- noise ratio in a given basis, 
and and not detected in another basis. For example, the wall is detected at more than GOcx by the 
ridgelet transform, and less than 5a by the wavelet transform. The line is detected almost at lOo" 
by the beamlet transform, and is under a 3a detection level using wavelets. These results shows 
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the importance of using several transforms for an optimal detection of all features contained in 
a data set. 

5.2 Experiment 2 
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Figure 14: Simulated data sets. Top, the Voronoi vertices point pattern (left) and the galaxies 
of the GIF A-CDM N-body simulation (right). The bottom panels show one 10 h^^ width slice 
of the each data set. 

We use here two simulated data sets to illustrate the discriminative power of the multiscale 
methods. The first one is a simulation from stochastic geometry. It is based on a Voronoi 
model. The second one is a mock catalog of the galaxy distribution drawn from a A-CDM N- 
body cosmological model[12_. Both processes have very similar two-point correlation functions 
at small scales, although they look quite different and have been generated following completely 
different algorithms. 

• the first comes from a Voronoi simulation: We locate a point in each of the vertices of a 
Voronoi tessellation of 1.500 cells defined by 1500 nuclei distributed following a binomial 
process. There are 10085 vertices lying within a box of 141.4 h~^ Mpc side. 

• the second point pattern represents the galaxy positions extracted from a cosmological 
A-CDM N-body simulation. The simulation has been carried out by the Virgo consortium 
and related groups (see |http : //www . mpa-gar ching . mpg . de/Virgo I . The simulation is a 
low-density (fi = 0.3) model with cosmological constant A = 0.7. It is, therefore, the 
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Figure 15: The two-point correlation function of the Voronoi vertices process and the GIF 
A-CDM N-body simulation. They are very similar in the range [0.02,2] h~^ Mpc. 

closer set to the real galaxy distribution^^. There are 15445 galaxies within a box with 
side 141.3 h~^ Mpc. Galaxies in this catalog have stellar masses exceeding 2 x 10^*^ Mq. 

Figure El shows the two simulated data set, and Figure El shows the two-point correlation 
function curve for the two point processes. The two point fields are different, but as it can be 
seen in Fig. 1151 both have very similar two-point correlation functions in a huge range of scales 
(2 decades). 

We have applied the three transforms to each data set, and we have calculated the skewness 



vector S 



i^S-un Sr, S"^ 



and the kurtosis vector K 



,,kl,kl) at each scale j. s- 









are 



respectively the skewness at scale j of the wavelet coefficients, the ridgelet coefficients and the 
beamlet coefficients. ki],kr, kl are respectively the kurtosis at scale j of the wavelet coefficients, 
the ridgelet coefficients and the beamlet coefficients. Figure El shows the kurtosis and the 
skewness vectors of the two data set at the two first scales. On the contrary to the two-point 
correlation function, this analysis shows strong differences between the two data set, particularly 
on the wavelet axis, which indicates that the second data contains more, or with a higher density, 
clusters than the first one. 

5.3 Experiment 3 

In this experiment, we have used a A-CDM simulation based on a N-body and hydrodynamical 
code, called RAMSES |3ZI- The code is based on Adaptive Mesh Refinement (AMR) technique, 
with a tree-based data structure allowing recursive grid refinements on a cell-by-cell basis. The 
simulated data have been obtained using 256^ particles and 4.1 x 10^ cells in the AMR grid, 
reaching a formal resolution of 8192^. The box size was set to lOO/i^l Mpc, with the following 
cosmological parameters: 



n. 



A 



0.3 n 

h = 0.7 as 



= 0.7 Qb = 0.039 
0.92 



(13) 
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Figure 16: Skewness and kurtosis for the two simulated data set. 

We used the results of this simulation at six different redshifts {z = 5, 3, 2, 1, 0.5, 0). Fig. El 
shows a projection of the simulated cubes along one axis. We have applied the 3D wavelet 
transform, the 3D beamlet transform and the 3D ridgelet transform on the six data set, and 
we calculate for each transform the standard deviation of the different scales. We will note 
'^Wz j^'^Rz j^'^B z j ^^^ variance of the scale j relative to the transformation of the cube at 
redshift z by respectively the wavelet, the ridgelet and the beamlet transform. 

Figure ITHl shows respectively from top to bottom the wavelet spectrum P^(z,j) = cfyy^ ■, 
the beamlet spectrum Pi,{z,j) = cr'^ ^ j and the ridgelet spectrum Pr{z,j) = crj^^j- In order 
to see the evolution of matter distribution with the redshift and scale, we calculate the ratio 
MUj,z) = ^ and MrUJ,z) = f^. 

Figure El shows the Mi and M2 curves as a function of z and Figure 1201 shows the Mi and 
M2 curves as a function of the scale number j. 

The Ml curve does not show too much evolution, while the M2 curve presents a significant 
slope. This shows that the beamlet transform is much more sensible to the formation of clusters 
than the ridgelet transform. This is not surprising since the beamlet function support is much 
smaller than the ridgelet function support. M2 increases with z showing clearly the cluster 
formation. This indicates that the combination of multiscale transformations allows us to get 
some information about the degree of clustering, filamentarity, and sheetedness. 
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Figure 17: A-CDM simulation at different redshifts. 

6 Conclusion 

We have introduced in this paper two new methods to analyze catalogs of galaxies. The first one 
consists in estimating the real underlying density though a wavelet denoising. Structures are 
first detected in the wavelet space, and an iterative reconstruction is performed. A smoothness 
constraint based on the li norm of the wavelet coefficients is used, which reduce the amount 
of artifact in the reconstructed density, especially the ringing artifacts around strong features 
which are due to the wavelet function shape. We have shown that such an approach leads to 
much more reliable results than a Gaussian filtering when we want to derive a Genus curve 
from the catalog. This could also be true for other techniques which require to pre-process the 
data with a Gaussian filtering. The wavelet denoising preserves the resolution of the detected 
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features whatever their sizes, and remove the noise in a non-hnear way at ah the scales. 

The second approach does not require to detect anything. It is based on the analyzing of the 
distribution of coefficients obtained by several niultiscale transforms. We have introduced two 
new multiscale decompositions, the 3D ridgelet transform and the 3D beamlet transform. We 
have described how to implement them using the FFT. Then we have shown that combining the 
information related to wavelet, ridgelet and beamlet coefficients leads to a new way to describe 
a data set. We have used in this paper the skewness and the kurtosis, but other recent statistic 
estimator such the Higher Criticism [^ could be used as well. Each multiscale transform is very 
sensible to one kind of feature, the wavelet to clusters, the beamlet to filaments, and the ridgelet 
to walls. A similar method has been proposed for analyzing CMB maps [22] where both the 
curvelet and the wavelet transform were used for the detection and the discrimination of non 
Gaussianities. This combined multiscale statistic is very powerful and we have shown that two 
data set that cannot be distinguished using a two point correlation function are clearly identified 
as different using our method. We believe that such an approach will permit to better constraint 
the cosmological models. 
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